1. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a punktami adresowymi.

maxpop = max(c(out_model_punkty_adresowe$TOT, out_model_punkty_adresowe$EST_POP))

p1 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = TOT)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) + 
  labs(title = "Ryc.1 Mapa oryginalnych wartości liczby ludności") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

p2 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = EST_POP)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) + 
  labs(title = "Ryc.2 Mapa estymowanych wartości liczby ludności") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p1)

print(p2)

p3 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = RES)) +
  scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(out_model_punkty_adresowe$RES), max(out_model_punkty_adresowe$RES))) + 
  labs(title = "Ryc.3 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p3)

out_model_punkty_adresowe$RES_reclass = case_when(
    out_model_punkty_adresowe$RES < 0 ~ "bledy ujemne",
    out_model_punkty_adresowe$RES == 0 ~ "brak bledu",
    out_model_punkty_adresowe$RES > 0 ~ "bledy dodatnie",
  )

p4 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = RES_reclass)) +
  scale_fill_manual(name = "Typ błędu", 
                    values = c("bledy dodatnie" = "red", "brak bledu" = "white", "bledy ujemne" = "blue")) +
  labs(title = "Ryc.4 Mapa zreklasyfikowanych błędów modelu") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p4)

meanerror = mean(out_model_punkty_adresowe$RES)
rmse = sqrt(mean((out_model_punkty_adresowe$RES)^2))

print(paste("średni błąd estymacji modelu wynosi", round(meanerror,2), "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.67 natomiast pierwiastek średniego błędu kwadratowego wynosi 30.44"

1.1. Ocena modelu

Maciejo,
YOU KNOW WHAT TO DO

2. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a liczbą mieszkań.

2.1 Liniowy model zależności między liczbą ludności a liczbą mieszkań w siatce 1km

# wykres rozrzutu
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

model_liniowy <- lm(TOT ~ N_BUDYNKI, data = pop1km)

par(mfrow = c(2, 2))
plot(model_liniowy, which = 1)
plot(model_liniowy, which = 2)
plot(model_liniowy, which = 3)
plot(model_liniowy, which = 4)
mtext("Ryc. 5 Wykresy diagnostyczne modelu liniowego", outer = TRUE, line = -1.5, cex = 1)

W jakim stopniu liczba ludności (TOT) jest wyjaśniana przez liczbę mieszkań?

2.2. Jakie są statystyki reszt? (residuals)

Ryc.6 Interaktywne wykresy reszt 
## Statystyki
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -98.6520  -2.5883  -0.7094   0.0000  -0.5883 174.7735
## Średnia:  3.738246e-16

3. Ocena dopasowania modelu (wykresy diagnostyczne, identyfikacja wartości odstających).

Wartości odstające można zidentyfikować wykorzystując wykresy diagnostyczne, a w szczególności wykres Residuals vs. Leverage.

p <- autoplot(model_liniowy) + theme_classic()

gridExtra::grid.arrange(grobs = p@plots, top = "Ryc. 7 Wykresy diagnostyczne modelu liniowego")

#dopasowanie modelu
smr_model_lm <-summary(model_liniowy)

c( R2 = smr_model_lm$r.squared, R2_adj = smr_model_lm$adj.r.squared)
##        R2    R2_adj 
## 0.7019419 0.6990481

~70% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI. 

# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy, level = 0.99)
##                 0.5 %   99.5 %
## (Intercept) -7.022728 8.441527
## N_BUDYNKI    3.275627 4.603257
outlierTest(model_liniowy)
##     rstudent unadjusted p-value Bonferroni p
## 21  9.742598         3.0248e-16   3.1760e-14
## 25  6.813338         6.7995e-10   7.1395e-08
## 34 -5.176442         1.1389e-06   1.1958e-04
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + 
  geom_point() +
  geom_point(data = pop1km[c(21, 25, 34),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
  labs(title = "Ryc.8 Wykres rozrzutu z oznaczonymi punktami o wartościach odstających") +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

meanerror = mean(model_liniowy$residuals)
rmse = sqrt(mean((model_liniowy$residuals)^2))

print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.73824648485021e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 26.25"

PYTANIA Czego dowiadujemy się z wykresów diagnostycznych? Czy model jest dobrze dopasowany? Czy spełnione są założenia regresji liniowej?

4. Model po usunięciu wartości odstających

4.1 Model regresji

# Usunięcie wierszy 21, 25 i 34
pop1km_n <- pop1km[-c(21, 25, 34), ]

# wykres rozrzutu
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
  labs(title = "Ryc.9 Wykres rozrzutu po usunięciu punktów z wartościami odstającymi") +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# model liniowy
model_liniowy_n <- lm(TOT ~ N_BUDYNKI, data = pop1km_n)

p1 <- autoplot(model_liniowy_n) + theme_classic()

gridExtra::grid.arrange(grobs = p1@plots, top = str_wrap("Ryc. 10 Wykresy diagnostyczne modelu liniowego po usunięciu punktów z wartościami odstającymi"))

4.2. Ocena modelu

#dopasowanie modelu
smr_model_lm_n <-summary(model_liniowy_n)

c( R2 = smr_model_lm_n$r.squared, R2_adj = smr_model_lm_n$adj.r.squared)
##        R2    R2_adj 
## 0.8914342 0.8903486

~89% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI. 

# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy_n, level = 0.99)
##                 0.5 %   99.5 %
## (Intercept) -2.413594 3.198192
## N_BUDYNKI    3.198037 3.843294
outlierTest(model_liniowy_n)
##     rstudent unadjusted p-value Bonferroni p
## 65 -6.324628         7.3790e-09   7.5266e-07
## 33  4.742867         7.0854e-06   7.2271e-04
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + 
  geom_point() +
  geom_point(data = pop1km_n[c(65, 33),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
  labs(title = str_wrap("Ryc.11 Wykres rozrzutu z oznaczonymi punktami o wartościach odstających (po usunięciu punktów)")) +
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'


#WIZUALIZACJA
par(mfrow = c(2, 4))
plot(model_liniowy)
plot(model_liniowy_n)
mtext(str_wrap("Ryc. 12 Wykresy diagnostyczne modelu liniowego (pierwszy wiersz - przed usunięciem, drugi wiersz - po)"), adj = 1, line = 20.5, font = 2, cex = 0.8)

meanerror = mean(model_liniowy_n$residuals)
rmse = sqrt(mean((model_liniowy_n$residuals)^2))

print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 2.30805383656037e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 9.14"

5. Wizualizacja wyników modelu na mapie

5.1. Mapa pokazująca oryginalne wartości liczby ludności (TOT)

maxpop1 = max(c(pop1km$TOT, pop1km$EST_POP))

ggplot() +
  geom_sf(data = pop1km, aes(fill = TOT)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop1)) +
  labs(fill = "Liczba ludności", title = "Ryc.13 Mapa oryginalnych wartości liczby ludności") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

5.2. Mapa pokazującą estymowaną liczbę ludności (EST_POP)

# Przewidywanie EST_POP na podstawie modelu liniowego
pop1km$EST_POP <- predict(model_liniowy, newdata = pop1km)

# Tworzenie mapy z estymowaną liczbą ludności (EST_POP)
ggplot() +
  geom_sf(data = pop1km, aes(fill = EST_POP)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop1)) +
  labs(fill = "Estymowana liczba ludności", title = "Ryc.14 Mapa estymowanych wartości liczby ludności") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

5.3. Mapa reszt (różnic między wartością obserwowaną TOT oraz EST_POP)

# Obliczenie reszt
pop1km$Reszty <- pop1km$TOT - pop1km$EST_POP

# Tworzenie mapy reszt
ggplot() +
  geom_sf(data = pop1km, aes(fill = Reszty)) +
  scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(pop1km$Reszty), max(pop1km$Reszty))) +
  labs(fill = "Reszty", title = "Ryc.15 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

5.4. Mapa reszt przeklasyfikowaną do 3 klas (błędy ujemne, brak błędu, błędy dodatnie)

# Przeklasyfikowanie reszt do trzech klas
pop1km$Klasy_reszt = case_when(
    pop1km$Reszty < 0 ~ "Niedoszacowane",
    pop1km$Reszty == 0 ~ "Brak błędu",
    pop1km$Reszty > 0 ~ "Przeszacowane",
  )

# Tworzenie mapy przeklasyfikowanych reszt
ggplot() +
  geom_sf(data = pop1km, aes(fill = Klasy_reszt)) +
  scale_fill_manual(name = "Typ błędu", values = c("Przeszacowane" = "red", "Brak błędu" = "white", "Niedoszacowane" = "blue")) +
  labs(title = "Ryc.16 Mapa zreklasyfikowanych błędów modelu") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

6. Statystyczna i przestrzenna analiza rokładu błędów modelu regresji na podstawie map oraz wyników modelu

## Statystyki reszt:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -98.6520  -2.5883  -0.7094   0.0000  -0.5883 174.7735
## Średnia reszt:  -5.399296e-16
## Odchylenie standardowe reszt:  26.37129
# Histogram reszt
ggplot(pop1km, aes(x=Reszty)) +
  geom_histogram(bins=30, fill="blue", color="black") +
  labs(title = "Ryc.17 Histogram reszt") +
  xlab("Reszty") +
  ylab("Częstość") +
  theme_minimal()

7. Mapa rozmieszczenia ludności w siatce 100m, dane pomocnicze - liczba mieszkań.

leaflet(data = pop_100m1, width = "100%", height = "600px") %>%
  addTiles() %>%
  addPolygons(
    data = pop_100m1,
    fillColor = ~pal(TOT),
    group = "Liczba ludności",
    color = "#BDBDC3",
    weight = 1,
    opacity = 1,
    fillOpacity = 1,
    highlightOptions = highlightOptions(
      weight = 2,
      color = "#666",
      fillOpacity = 1,
      bringToFront = TRUE
    ),
    label = ~paste("Liczba ludności:", TOT),
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = pal,
    values = ~TOT,
    opacity = 1,
    title = "Liczba ludności",
    position = "bottomright"
  ) %>%
  addLayersControl(
    baseGroups = c("Podkład mapy"),
    overlayGroups = c("Liczba ludności"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  setView(lng = mean(st_coordinates(pop_100m1)[,1]), 
          lat = mean(st_coordinates(pop_100m1)[,2]), 
          zoom = 12)

8. Porównanie wyników obu modeli